#include<fstream>
#include<cmath>
#include<algorithm>
#include"system_k.h"
#include"help.h"

using std::ifstream;
using std::abs;
using std::swap;

/************************************************************/
//!!!!TODO Orbin>>complex.real() is not c++ standard!!!!!
// void System_Tz_WS::setupOrbitals()
// {
//   ifstream OrbIn( OrbFile.c_str() );
//   if(!OrbIn) std::cerr<<"can't open Orbital Input file!\n";
//   int totalOrbitals;
//   OrbIn.ignore(128,'\n');
//   OrbIn.ignore(128,'\n');
//   OrbIn.ignore(128,'\n');
//   OrbIn.ignore(128,'\n');
//   OrbIn.ignore(128,'\n');
//   OrbIn.ignore(128,'\n');
//   OrbIn.ignore(128,'\n');
//   OrbIn.ignore(128,'\n');
//   OrbIn.ignore(45,'\n');
//   OrbIn>>totalOrbitals;
//   OrbIn.ignore(128,'\n');
//   OrbIn.ignore(128,'\n');
//   Orbitals.resize(totalOrbitals);
//   int num,n,l,jj,tz,nnl;
//   double e_real,e_imag,eval_real,eval_imag;
//   for(int i=0;i<totalOrbitals;i++)
//     {
//       OrbIn.ignore(128,':');
//       OrbIn>>num>>n>>l>>jj>>tz>>nnl>>e_real>>e_imag>>eval_real>>eval_imag;
//       OrbIn.ignore(128,'\n');
//       Orbitals[i]=OrbitalType(l,jj,tz,complexd(e_real,e_imag));
//     }
//   setupGroups();
//   OneBodyMat.resize( totalOrbitals,vector<complexd>(totalOrbitals) );
//   for(int i=0;i<totalOrbitals;i++)
//     for(int j=0;j<totalOrbitals;j++)
//       {
// 	double real,imag;
// 	OrbIn>>real>>imag;
// 	OneBodyMat[i][j]=complexd(real,imag);
//       }
// }

// const int System_Tz_WS::MatEledim=1;//initialize of static variable
// void System_Tz_WS::setupTwoBodyMat()
// {
//   setupIndiceIn2BStates();
//   int totalChannels=TwoBodyChannels.size();
//   TwoBodyMat.resize(totalChannels);
//   for(int channel=0;channel<totalChannels;channel++)
//     {
//       int num2BInChannel=TwoBodyChannels[channel].size();
//       TwoBodyMat[channel].clear();
//       TwoBodyMat[channel].resize( (num2BInChannel*(num2BInChannel+1))/2 );
//     }
//   ifstream MatIn(MatFile.c_str());
//   if(!MatIn) std::cerr<<"can't open Mat Input file!\n";
//   MatIn.ignore(128,'\n');
//   MatIn.ignore(128,'\n');
//   MatIn.ignore(128,'\n');
//   MatIn.ignore(128,'\n');
//   MatIn.ignore(128,'\n');
//   MatIn.ignore(128,':');

//   int totalLines,ppLines,pnLines,nnLines;
//   MatIn>>totalLines>>ppLines>>pnLines>>nnLines;
//   MatIn.ignore(128,'\n');
//   MatIn.ignore(128,'\n');
//   MatIn.ignore(128,'\n');
//   int Tz,Par,JJ,a,b,c,d;
//   TwoBodyMatEle MatEle(MatEledim);
//   for(int i=0;i<ppLines;i++)
//     {
//       MatIn>>Tz>>Par>>JJ>>a>>b>>c>>d;
//       for(int j=0;j<MatEledim;j++)
// 	MatIn>>MatEle[j].real()>>MatEle[j].imag();
//       double t1;
//       MatIn>>t1>>t1>>t1>>t1>>t1>>t1;
//       --a;--b;--c;--d;//the indice start from zero here, the input interaction from jensen start from 1.
//       set2BmatAt(a,b,c,d,JJ/2,MatEle);
//     }
//   for(int i=0;i<pnLines;i++)
//     {
//       MatIn>>Tz>>Par>>JJ>>a>>b>>c>>d;
//       for(int j=0;j<MatEledim;j++)
// 	MatIn>>MatEle[j].real()>>MatEle[j].imag();
//       double t1;
//       MatIn>>t1>>t1>>t1>>t1>>t1>>t1;
//       int phase=1;
//       --a;--b;--c;--d;
//       if(a>b)
// 	{
// 	  swap(a,b);
// 	  phase*=( (Orbitals[a].jj+Orbitals[b].jj-JJ)/2 + 1 )%2?-1:1;
// 	}
//       if(c>d)
// 	{
// 	  swap(c,d);
// 	  phase*=( (Orbitals[c].jj+Orbitals[d].jj-JJ)/2 + 1 )%2?-1:1;
// 	}
//       MatEle*=phase;
//       set2BmatAt(a,b,c,d,JJ/2,MatEle);
//     }
//   for(int i=0;i<nnLines;i++)
//     {
//       MatIn>>Tz>>Par>>JJ>>a>>b>>c>>d;
//       for(int j=0;j<MatEledim;j++)
// 	MatIn>>MatEle[j].real()>>MatEle[j].imag();
//       double t1;
//       MatIn>>t1>>t1>>t1>>t1>>t1>>t1;
//       --a;--b;--c;--d;
//       set2BmatAt(a,b,c,d,JJ/2,MatEle);
//     }
// }

// complexd System_Tz_WS::get1B(int A,int bra,int ket) const
// {
//   complexd onebody=OneBodyMat[bra][ket];
//   double fac=(1-1./A)*hbar_c*hbar_c*0.5/mc2;
//   onebody*=fac;
//   return onebody;
// }

// complexd System_Tz_WS::get2B(int A,int a,int b,int c,int d,int J) const
// {
//   TwoBodyMatEle MatEle(System_Tz_WS::MatEledim);
//   get2BmatAt(a,b,c,d,J,MatEle);
//   if(MatEle.empty()) return 0;
//   complexd V=MatEle[0];
//   const double sqrt2=sqrt(2.);
//   V *= ( (a==b)? sqrt2 : 1.0 );
//   V *= ( (c==d)? sqrt2 : 1.0 );
//   return V;
// }

//****************************************************************************************

void System_Tz_K::setupOrbitals()
{
  ifstream OrbIn("./sppn.dat");
  if(!OrbIn) std::cerr<<"can't open Orbital Input file!\n";

  double lab_cutoff=10;
  gauleg(0,lab_cutoff,labk,labwk);

  // gauleg(0,1,labk,labwk);
  // for(int i=0;i<labk.size();i++)
  //   {
  //     labwk[i]=labwk[i]/(labk[i]*labk[i]);
  //     labk[i]=1/labk[i]-1;
  //   }
  int totalOrbitals;
  OrbIn.ignore(128,'\n');
  OrbIn.ignore(128,'\n');
  OrbIn.ignore(128,'\n');
  OrbIn.ignore(128,'\n');
  OrbIn.ignore(128,'\n');
  OrbIn.ignore(128,'\n');
  OrbIn.ignore(128,'\n');
  OrbIn.ignore(128,'\n');
  OrbIn.ignore(45,':');
  OrbIn>>totalOrbitals;
  OrbIn.ignore(128,'\n');
  OrbIn.ignore(128,'\n');
  Orbitals.resize(totalOrbitals*meshPoints);
  int num,n,l,jj,tz,nnl;
  double e,eval;
  double k,wk;
  for(int i=0;i<totalOrbitals;i++)
    {
      OrbIn.ignore(128,':');
      OrbIn>>num>>n>>l>>jj>>tz>>nnl;
      OrbIn.ignore(128,'\n');
      for(int j=0;j<meshPoints;j++)
	{
	  k=labk[j];
	  wk=labwk[j];
	  e=k*k*0.5*hbar_c*hbar_c/mc2;
	  Orbitals[i*meshPoints+j]=OrbitalType(l,jj,tz,e,k,wk);
	}
    }
  setupGroups();
}

void System_Tz_K::setupTwoBodyMat()
{
}

double System_Tz_K::get1B(int A,int bra,int ket) const
{
  double t=Orbitals[bra].k*Orbitals[bra].k*0.5*hbar_c*hbar_c/mc2;
  return (bra==ket)?(1-1./A)*t:0;
}

double System_Tz_K::get2B(int a,int b,int c,int d,int J) const
{	  
  int la=Orbitals[a].l;
  int jja=Orbitals[a].jj;
  double pa=Orbitals[a].k;
  double wa=Orbitals[a].wk;

  int lb=Orbitals[b].l;
  int jjb=Orbitals[b].jj;
  double pb=Orbitals[b].k;
  double wb=Orbitals[b].wk;

  int lc=Orbitals[c].l;
  int jjc=Orbitals[c].jj;
  double pc=Orbitals[c].k;
  double wc=Orbitals[c].wk;
	  
  int ld=Orbitals[d].l;
  int jjd=Orbitals[d].jj;
  double pd=Orbitals[d].k;
  double wd=Orbitals[d].wk;
	  
  int isospin_tzab=Orbitals[a].tz+Orbitals[b].tz;
  int isospin_tzcd=Orbitals[c].tz+Orbitals[d].tz;
  if(isospin_tzab!=isospin_tzcd) return 0;
  if((la+lb+lc+ld)%2) return 0;
  if(J> min(jja+jjb,jjc+jjd)/2) return 0;
  if(J< max(abs(jja-jjb),abs(jjc-jjd))/2) return 0;
  int pha=1;
  if(Orbitals[a].tz>Orbitals[b].tz)
    {
      swap(la,lb);
      swap(jja,jjb);
      swap(pa,pb);
      pha*=phase( (jja+jjb)/2-J+1 );
    }
  if(Orbitals[c].tz>Orbitals[d].tz)
    {
      swap(lc,ld);
      swap(jjc,jjd);
      swap(pc,pd);
      pha*=phase( (jjc+jjd)/2-J+1 );
    }
  double t=v_k(la,jja,pa,lb,jjb,pb,lc,jjc,pc,ld,jjd,pd,isospin_tzab/2,J)*sqrt(wa*wb*wc*wd);
  return t*pha;
}

void System_Tz_K::printOrbitals()
{
  for(int i=0;i<Orbitals.size();i++)
    {
      std::cout<<Orbitals[i].l<<"\t"<<Orbitals[i].jj<<"\t"<<Orbitals[i].tz<<"\t"<<Orbitals[i].k<<"\t"<<Orbitals[i].wk<<std::endl;
    }
}
